;; 批量处理CALIPSO VFM数据并生成地理分布图
; 支持处理指定目录下的所有.hdf文件

; 输入目录（包含HDF文件） G:\CALIPSO\code\python\su\calipso-vfm-visualization-main\data
input_dir = 'G:\CALIPSO\code\python\su\calipso-vfm-visualization-main\data\'

; 输出目录   G:\CALIPSO\PRO_CODE\output_geo_wrold_track\
output_dir = 'G:\CALIPSO\PRO_CODE\output_geo_wrold_track\\'

; 确保输出目录存在
IF ~FILE_TEST(output_dir, /DIRECTORY) THEN FILE_MKDIR, output_dir

; 查找所有.hdf文件
file_pattern = input_dir + '*.hdf'
file_list = FILE_SEARCH(file_pattern)
n_files = N_ELEMENTS(file_list)

; 检查是否找到文件
IF n_files EQ 0 THEN BEGIN
  PRINT, '未找到HDF文件！请检查输入目录：', input_dir
  RETURN
ENDIF

PRINT, '找到', n_files, '个HDF文件，开始处理...'

; 批量处理每个文件
FOR i = 0, n_files-1 DO BEGIN
  file_name = file_list[i]
  PRINT, '处理文件:', file_basename(file_name)
  base_name = file_basename(file_name)
  file_core = strmid(base_name, 0, strlen(base_name) - 4)
  ; 打开文件
  fid = HDF_SD_START(file_name, /READ)
  IF fid LT 0 THEN BEGIN
    PRINT, '错误：无法打开文件', file_name
    CONTINUE
  ENDIF

  ; Define data name. Use "HDFView" or "hdp" tool to check datasets
  ; available inside the file.
  data_name = "Feature_Classification_Flags"
  index = HDF_SD_NAMETOINDEX(fid, data_name)

  ; Retrieve FCF data.
  sds_id = HDF_SD_SELECT(fid, index)
  HDF_SD_GETDATA, sds_id, fcf

  ; Retrieve lat data.
  data_name = "Latitude"
  index = HDF_SD_NAMETOINDEX(fid, data_name)
  sds_id = HDF_SD_SELECT(fid, index)
  HDF_SD_GETDATA, sds_id, lat

  ; Retrieve lon data.
  data_name = "Longitude"
  index = HDF_SD_NAMETOINDEX(fid, data_name)
  sds_id = HDF_SD_SELECT(fid, index)
  HDF_SD_GETDATA, sds_id, lon

  ; End access to SDS.
  HDF_SD_ENDACCESS, sds_id

  ; Close file.
  HDF_SD_END, fid

  ; Users need to understand the layout of the Feature_Classification_Flag
  ; dataset.
  ; ... 省略重复注释 ...

  ; Subset data at 2500m (= -0.5km + 30m * 100) altitude.
  alt_index = 1256
  data = fcf(alt_index, *)

  ; Select the first 3 bit for the feature type.
  data = data AND '0007'X

  ; lat/lon are 2-D, Nx1 array.
  lat = reform(lat)
  lon = reform(lon)
  data = reform(data)
  dim = size(data, /dim)
  longname = 'Feature Type at altitude=2500m'

  ; Define the color table.
  levels = 8
ct = COLORTABLE([[150, 0,   255,  0,   255,  200, 100, 50], $
                 [150, 0,   255,  255, 0,    100, 50,  25], $
                 [150, 255, 0,    0,   0,    255, 255, 125]], $
                 NCOLORS = levels, /TRANSPOSE)
;  ct = COLORTABLE([[255, 224, 135, 255, 255, 46,  152, 0], $
;    [255, 255, 206, 215,  0,  139, 251, 0], $
;    [255, 255, 235,  0,   0,  87,  152, 0]], $
;    NCOLORS = levels,   /TRANSPOSE)

  ; Generate the plot.
  m = MAP('Geographic', TITLE=file_core, FONT_SIZE=9, /BUFFER)
  t1 = TEXT(0.35, 0.2, longname)

  ; MAGNITUDE requires BYTE type data.
  data = BYTE(data)

  ; We use SCATTERPLOT because data is 2-d lat/lon swath.
  c1 = SCATTERPLOT(lon, lat, OVERPLOT=m, $
                    MAGNITUDE=data, $
                    RGB_TABLE=ct, $
                    POSITION=[0.1, 0.1, 0.83, 0.9],$
                    /SYM_FILLED, SYMBOL='o', SYM_SIZE=0.1)
                  mc = MAPCONTINENTS()

  ; We need a custom colorbar because we use SCATTERPLOT().
  cb = COLORBAR(RGB_TABLE=ct, RANGE=[0,8], ORIENTATION=1, BORDER=1,$
                TICKVALUES=[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5], $
                TICKNAME=['invalid', 'clear', 'cloud', 'aerosol', $
                'strato.', 'surface', 'subsurf.', 'no signal'], $
                TEXTPOS=1, POSITION=[0.85,0.2,0.87,0.8], TITLE=units)

  ; 生成输出文件名
    output_file = output_dir + 'geo_wrold_track_' + file_core + '.png'

  ; 保存图像
  ; 图像尺寸（英寸）
  image_width = 10.0
  image_height = 8.0

  ; 固定输出分辨率（DPI）
  output_dpi = 600
  
  ; 计算像素尺寸（根据固定的600 DPI）
  pixel_width = ROUND(image_width * output_dpi)
  pixel_height = ROUND(image_height * output_dpi)
  
  ; 保存图像并设置像素尺寸（固定600 DPI）
  c1.save, output_file, HEIGHT=pixel_height, WIDTH=pixel_width
  PRINT, '已保存:', output_file, ' (', pixel_width, 'x', pixel_height, '像素, 600 DPI)'
  
;  c1.save, output_file, HEIGHT=600, WIDTH=800
;  ; 保存图像并设置像素尺寸（控制分辨率）
;  c1.save, output_file, HEIGHT=pixel_height, WIDTH=pixel_width
;  PRINT, '已保存:', output_file, ' (', pixel_width, 'x', pixel_height, '像素,', output_dpi, 'DPI)'
;  PRINT, '已保存:', output_file
ENDFOR

PRINT, '批量处理完成！'
end